GRCh38 vs hg38

• GRCh38: The Human Reference Genome maintained by the Genome Reference Consortium (GRC) is the 38th edition of the Genome Reference Consortium Human series.

• hg38: indicates the GRCh38 version used by UCSC

• Same

  GRCh38 and hg38 share the same genome sequence, so most of the coordinates are identical

• Different

  UCSC uses the UCSC-specific mitochondrial genome (chrM) in hg38, while GRCh38 may use a different mitochondrial version (such as rCRS or GRCh38_MT).

  UCSC has added alternate haplotypes in hg38

What is the Haplotype Patches or alternate haplotypes ?

• 单倍型补丁（Haplotype Patches）

  chr1    1000000 1000500 region_A    100 +
  chr1    2000000 2000600 region_B    200 -
  chr1_KI270706v1_random  50000   50500   patch_region    150 +

  • chr1_KI270706v1_random 是 单倍型补丁，用于修复 chr1 上某个复杂区域的问题。

  • 这些补丁的坐标通常是独立的，它们不会直接与主染色体坐标系对应（即 chr1 和 chr1_KI270706v1_random 可能是相互独立的）

• 替代单倍型（Alternate Haplotypes）

  chr6    3000000 3000500 MHC_region_ref  500 +
  chr6_GL000256v2_alt 3000000 3000500 MHC_region_alt  500 +

  • chr6 是标准染色体的 MHC （Major Histocompatibility Complex，主要组织相容性复合体）(负责编码细胞表面的分子，用于抗原呈递和免疫系统识别)（主参考序列）。

  • chr6_GL000256v2_alt 是 MHC 替代单倍型（alt haplotype），用于提供该区域的另一种版本，表示在不同人群中可能存在的变异。

mRatBN7.2 vs rn7

It’s the same version of the genome, just from a different database

• mRatBN7.2 from Ensemble

• rn7 from UCSC

Conclusion

GRanges objects

• Each GRanges object can include the following

  • components:Chromosome Names (seqnames)

  • Start and End Positions (ranges)

  • Strand Information

    • Positive strand: +

    • Negative strand: -

    • Unspecified strand: *

  • Additional Metadata(score, Annotations,Other custom attributes)

# loading package
# processing genomic data.
library(rtracklayer)

## 载入需要的程序包：GenomicRanges

## 载入需要的程序包：stats4

## 载入需要的程序包：BiocGenerics

## 
## 载入程序包：'BiocGenerics'

## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs

## The following objects are masked from 'package:base':
## 
##     anyDuplicated, aperm, append, as.data.frame, basename, cbind,
##     colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
##     get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
##     match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
##     Position, rank, rbind, Reduce, rownames, sapply, saveRDS, setdiff,
##     table, tapply, union, unique, unsplit, which.max, which.min

## 载入需要的程序包：S4Vectors

## 
## 载入程序包：'S4Vectors'

## The following object is masked from 'package:utils':
## 
##     findMatches

## The following objects are masked from 'package:base':
## 
##     expand.grid, I, unname

## 载入需要的程序包：IRanges

## 
## 载入程序包：'IRanges'

## The following object is masked from 'package:grDevices':
## 
##     windows

## 载入需要的程序包：GenomeInfoDb

# crea two GRanges object
query <- GRanges(seqnames = c("chr1","chr1","chr1"),
                  ranges = IRanges(start = c(100, 300, 500), 
                                   end = c(200, 400, 600),
                                   names=c("a","b","c")),
                  score = c(10, 20, 30),
                  strand="+")

subject <- GRanges(seqnames = c("chr1", "chr1", "chr2"),
                  ranges = IRanges(start = c(150, 175, 550), 
                                   end = c(250, 450, 650),
                                   names=c("A","B","C")),
                  score = c(20, 40, 20),
                  strand="+")

# check the GRanges object
query

## GRanges object with 3 ranges and 1 metadata column:
##     seqnames    ranges strand |     score
##        <Rle> <IRanges>  <Rle> | <numeric>
##   a     chr1   100-200      + |        10
##   b     chr1   300-400      + |        20
##   c     chr1   500-600      + |        30
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

subject

## GRanges object with 3 ranges and 1 metadata column:
##     seqnames    ranges strand |     score
##        <Rle> <IRanges>  <Rle> | <numeric>
##   A     chr1   150-250      + |        20
##   B     chr1   175-450      + |        40
##   C     chr2   550-650      + |        20
##   -------
##   seqinfo: 2 sequences from an unspecified genome; no seqlengths

• You can use bed file as input to creat GRanges objects

IRanges objects

• IRanges(IntegerRanges) lack the informations of seqnames and strand and add the width informations compareing to GRanges

  # 创建 IRanges 对象
  ir <- IRanges(start = c(1, 3, 20), end = c(13, 7, 25))
  
  # 查看对象
  ir

  ## IRanges object with 3 ranges and 0 metadata columns:
  ##           start       end     width
  ##       <integer> <integer> <integer>
  ##   [1]         1        13        13
  ##   [2]         3         7         5
  ##   [3]        20        25         6

TxDb object

• TxDb (Transcript Database) is a object provided by Bioconductor for storing annotation information of genomes

• you can download from Bioconductor

  # loading packages
  library(AnnotationHub)

  ## 载入需要的程序包：BiocFileCache

  ## 载入需要的程序包：dbplyr

  ## 
  ## 载入程序包：'AnnotationHub'

  ## The following object is masked from 'package:rtracklayer':
  ## 
  ##     hubUrl

  # downland AnnotationHub database
  AH <- AnnotationHub()
  
  # Query all available Txdb databases
  query(AH , "Txdb.Rnorvegicus.UCSC.rn6")

  ## AnnotationHub with 6 records
  ## # snapshotDate(): 2024-10-28
  ## # $dataprovider: UCSC
  ## # $species: Rattus norvegicus
  ## # $rdataclass: TxDb
  ## # additional mcols(): taxonomyid, genome, description,
  ## #   coordinate_1_based, maintainer, rdatadateadded, preparerclass, tags,
  ## #   rdatapath, sourceurl, sourcetype 
  ## # retrieve records with, e.g., 'object[["AH52270"]]' 
  ## 
  ##             title                                      
  ##   AH52270 | TxDb.Rnorvegicus.UCSC.rn6.refGene.sqlite   
  ##   AH53781 | TxDb.Rnorvegicus.UCSC.rn6.refGene.sqlite   
  ##   AH61799 | TxDb.Rnorvegicus.UCSC.rn6.refGene.sqlite   
  ##   AH66179 | TxDb.Rnorvegicus.UCSC.rn6.refGene.sqlite   
  ##   AH70598 | TxDb.Rnorvegicus.UCSC.rn6.refGene.sqlite   
  ##   AH84143 | TxDb.Rnorvegicus.UCSC.rn6.ncbiRefSeq.sqlite

  # According to AH ID to downland database
  # A larger AH ID indicates a newer version.
  TxDb <- AH[["AH70598"]]

  ## loading from cache

  ## 载入需要的程序包：GenomicFeatures

  ## 载入需要的程序包：AnnotationDbi

  ## 载入需要的程序包：Biobase

  ## Welcome to Bioconductor
  ## 
  ##     Vignettes contain introductory material; view with
  ##     'browseVignettes()'. To cite Bioconductor, see
  ##     'citation("Biobase")', and for packages 'citation("pkgname")'.

  ## 
  ## 载入程序包：'Biobase'

  ## The following object is masked from 'package:AnnotationHub':
  ## 
  ##     cache

• you can use GFF/GTF to build TxDb object

  

EnsDb

• 我们还可以使用 Ensembl-based EnsDb 注释数据库（由 ensembldb 包提供）作为TxDb objects 进行输入

• 但需要用seqlevelsStyle() 对染色体名称进行转换

  特性	EnsDb	TxDb
  数据来源	Ensembl数据库的基因组注释数据	多种注释源（如UCSC、NCBI）
  数据内容	基因、转录本、外显子、CDS等全面信息	主要包含转录本信息，如外显子、内含子
  适用场景	基因组级别的分析，特别是Ensembl注释	转录本分析，关注转录本结构与功能

GFF(GFF3)

• GFF(General Feature Format)是一种用来描述基因组特征的文件

• 现在我们所使用的大部分都是第三版（GFF3）

  ctg123 . mRNA            1300  9000  .  +  .  ID=mrna0001;Name=sonichedgehog
  ctg123 . exon            1300  1500  .  +  .  ID=exon00001;Parent=mrna0001
  ctg123 . exon            1050  1500  .  +  .  ID=exon00002;Parent=mrna0001
  ctg123 . exon            3000  3902  .  +  .  ID=exon00003;Parent=mrna0001
  ctg123 . exon            5000  5500  .  +  .  ID=exon00004;Parent=mrna0001
  ctg123 . exon            7000  9000  .  +  .  ID=exon00005;Parent=mrna0001

  列名	描述
  seqid	参考序列的 ID
  source	该特征的来源，未知时用 . 代替
  type	特征类型，如 gene、exon、CDS 等，建议使用 SO 术语
  start	起始坐标（1-based）
  end	终止坐标
  score	量化得分，空值用 . 代替
  strand	+ 表示正链，- 表示负链，. 表示无方向要求
  phase	CDS 需要此列，当前 CDS 片段的起始碱基在密码子中的偏移量。它的值可以是 0、1 或 2
  attributes	属性列表，格式为 key=value，不同属性用 ; 分隔

GTF (GTF2)

• GTF(gene transfer format)，主要是用来对基因进行注释

• 当前所广泛使用的GTF格式为第二版（GTF2）

  AB000381 Twinscan  CDS          700   707   .   +   2  gene_id "AB000381.000"; transcript_id "AB000381.000.1";
  AB000381 Twinscan  exon         900  1000   .   +   .  gene_id "AB000381.000"; transcript_id "AB000381.000.1";
  AB000381 Twinscan  start_codon  380   382   .   +   0  gene_id "AB000381.000"; transcript_id "AB000381.000.1";
  AB000381 Twinscan  stop_codon   708   710   .   +   0  gene_id "AB000381.000"; transcript_id "AB000381.000.1";

  列名	描述
  seqname	序列 ID，通常是染色体 ID 或 Contig ID
  source	该特征的来源，未知时用 . 代替
  type	特征类型，特征类型，如 CDS、start_codon(起始密码子)、stop_codon（终止密码子） （这一列必须注明）
  start	起始坐标（1-based）
  end	终止坐标
  score	量化得分，空值用 . 代替
  strand	+ 表示正链，- 表示负链，. 表示无方向要求
  frame	读码框（0、1 或 2，仅适用于 CDS）
  attributes	其他信息（键值对，key "value"; 格式）

Different from GFF of GTF

Loading chain file

• UCSC chain file download web (Recommend)

  (choose species ——choose genome version—–choose LifeOver files)

• Ensembl chain file download web (Only have human and mouse)

# decompressing chain file
if (file.exists("rn6ToRn7.over.chain")) {
  print("文件解压成功！")
} else {
  gunzip("rn6ToRn7.over.chain.gz")
  print("文件重新解压完成。")
}

## [1] "文件解压成功！"

# read chain file
chain_file <- import.chain("rn6ToRn7.over.chain")

Header Line

chain score tName tSize tStrand tStart tEnd qName qSize qStrand qStart qEnd id 

• score – chain score （How well the two sequences allign）

• tName – chromosome (reference/target sequence)

• tSize – chromosome size (reference/target sequence)

• tStrand – strand (reference/target sequence)

• tStart – alignment start position (reference/target sequence)

• tEnd – alignment end position (reference/target sequence)

• qName – chromosome (query sequence)

• qSize – chromosome size (query sequence)

• qStrand – strand (query sequence)

• qStart – alignment start position (query sequence)

• qEnd – alignment end position (query sequence)

• id – chain ID

Alignment Data Lines

size dt dq 

• Alignment Data Lines 将 chain分成许多个 block(一致的区域)(对齐块)，每一行代表一个block

• size – the size of the ungapped alignment (block)

• dt – the difference between the end of this block and the beginning of the next block (reference/target sequence)

• dq – the difference between the end of this block and the beginning of the next block (query sequence)

     chain 4900 chrY 58368225 + 25985406 25985566 chr5 151006098 - 43549808 43549970 2
     16      0       2
     60      4       0
     10      0       4
     70 

Per setting

# setting distance_threshold (base)
distance_threshold <- 30

# setting minimum width of peaks (base)
minimum_width <- 30

Function of merge peaks by group

merge_peaks_by_group <- function(df, distance_threshold) {
  df %>%
    group_by(group) %>%  # 按 group 分组
    filter(n_distinct(seqnames) == 1) %>%  # 只保留同一组内染色体相同的行
    arrange(group, start) %>%  # 按照 start 位置排序
    mutate(
      dist_to_prev = start - lag(end, default = first(end))  # 增加dist_to_prev这列，其为当前 start 和前一个 peak 的 end 之间的距离
    ) %>%
    ungroup() %>%  # 取消分组，以便下一步操作
    mutate(
      merge_group = cumsum(dist_to_prev > distance_threshold | is.na(dist_to_prev))  # 增加merge_group列，如果距离大于阈值，增加1（不合并）
    ) %>%
    group_by(group, merge_group) %>%  # 按照 group 和新的merge_group进行分组
    summarise(
      group_name = first(group_name),  # 保留第一个 group_name
      seqnames = first(seqnames),
      start = min(start),  # 取最小的 start
      end = max(end),  # 取最大的 end
      width = end - start,  # 计算新的 width
      strand = first(strand),  # 保留第一个 strand
      score = first(score),  # 保留第一个 score
      .groups = "drop"  # 去掉分组结构
    )
}

merge_peaks_subject <- merge_peaks_by_group(subject_data,distance_threshold)

Remove narrow width of peaks

• Evaluate the change in width after liftOver to assess the quality of merged peaks

• peak whose width increases too much after lifover is removed

# Extract original width
before_width <- bed_data[,c(22,4)]
colnames(before_width) <- c("group_name","before_width")


# merge original width 
merge_peaks_subject <- merge(merge_peaks_subject,before_width,by.x="group_name",by.y="group_name")


# calculate the change in width after liftOver to assess the quality
merge_peaks_subject$width_change <- merge_peaks_subject$width/merge_peaks_subject$before_width*100


# remove narrow width of peaks
merge_peaks_subject <- filter(merge_peaks_subject,merge_peaks_subject$width>=minimum_width)

findOverlaps()

• We use findOverlaps() to find the overlaps

# creat two GRanges object
query <- GRanges(seqnames = c("chr1","chr1","chr1"),
                  ranges = IRanges(start = c(100, 300, 500), 
                                   end = c(200, 400, 600),
                                   names=c("a","b","c")),
                  score = c(10, 20, 30),
                  strand="+")

subject <- GRanges(seqnames = c("chr1", "chr1", "chr2"),
                  ranges = IRanges(start = c(150, 175, 550), 
                                   end = c(250, 450, 650),
                                   names=c("A","B","C")),
                  score = c(20, 40, 20),
                  strand="+")

hits <- findOverlaps(query, #  query GRanges object
                         subject, #  subject GRanges object
                         type="any",
                         #minoverlap=1, # 两个区域之间的最小重叠碱基数
                         #maxgap=0, # 允许重叠区域之间的最大间隙
                         ignore.strand= F, # F means strand ising considered
                         select="all",
                         # invert=TRUE #If True, keep only the ranges in x that do not overlap ranges
                         )

• type

  • any: 认为query和subject存在任意的重合即认为发生了重叠（deafult）
  • within: 认为query只有完全在subject内才认为是重叠
  • start: 认为query中的start位置在subjectt内才认为是重叠
  • end: 认为query中的end位置在subject内才认为是重叠
  • equal：认为query中的start，end与subject完全相同才认为是重叠
  • When type is “any”, at least one of maxgap and minoverlap must be set to its default value.

  

• maxgap and minoverlap

  • Positive Values of maxgap：if maxgap = 50, two ranges can be up to 50 base pairs apart and still be considered overlapping

  • Negative Values of maxgap： if maxgap = -20, the two ranges must overlap by at least 20 base pairs.

  • minoverlap：Only ranges with a minimum of minoverlap overlapping positions are considered to be overlapping.

  

• select

  • all: Return all overlaps (default).

  • first: Return the first overlap for each query range

  • last: Return the last overlap for each query range

  • arbitrary: Return one arbitrary overlap for each query range

hits objects

• The function of findoverlaps() function will creat hits objects

hits

## Hits object with 3 hits and 0 metadata columns:
##       queryHits subjectHits
##       <integer>   <integer>
##   [1]         1           1
##   [2]         1           2
##   [3]         2           2
##   -------
##   queryLength: 3 / subjectLength: 3

• We can use countQueryHits()extracting numbers of overlap of query range

• Using countSubjectHits() extracting numbers of overlap of subject range

# Query
countQueryHits(hits)

## [1] 2 1 0

setNames(countQueryHits(hits), seqnames(subject))

## chr1 chr1 chr2 
##    2    1    0

# subject
countSubjectHits(hits)

## [1] 1 2 0

setNames(countSubjectHits(hits), seqnames(subject))

## chr1 chr1 chr2 
##    1    2    0

• Using overlapsRanges() to check the Specific range of overlaps

GenomicRanges::intersect(query, subject)

## GRanges object with 2 ranges and 0 metadata columns:
##       seqnames    ranges strand
##          <Rle> <IRanges>  <Rle>
##   [1]     chr1   150-200      +
##   [2]     chr1   300-400      +
##   -------
##   seqinfo: 2 sequences from an unspecified genome; no seqlengths

countOverlaps()

• countOverlaps()可以直接查看query与subject重叠的数目

• countOverlaps() 占用的内存非常少，再快速查看重叠数目非常有用

• countOverlaps() 用法与findOverlaps 一样

overlap_counts <- countOverlaps(query, subject)
overlap_counts

## a b c 
## 2 1 0

Extract the overlapping peaks

• The function of pintersect() find the overlapping portion of each range in one set with its paired range in other set

# Extract overlapping ranges
overlap_ranges <- pintersect(query[queryHits(hits)], 
                             subject[subjectHits(hits)])

# caulating the mean score
mean_scores <- rowMeans(cbind(
  query$score[queryHits(hits)],
  subject$score[subjectHits(hits)]
))

# # Add the mean scores as metadata to the overlap ranges
overlap_ranges$score <- mean_scores

overlap_ranges

## GRanges object with 3 ranges and 2 metadata columns:
##     seqnames    ranges strand |     score       hit
##        <Rle> <IRanges>  <Rle> | <numeric> <logical>
##   a     chr1   150-200      + |        15      TRUE
##   a     chr1   175-200      + |        25      TRUE
##   b     chr1   300-400      + |        30      TRUE
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

Find Overlaps with ≥50%(threshold)

# setting the threshold
threshold=0.5


# exact the overlap_widths
overlap_widths <- width(overlap_ranges)


# Calculate widths of query and subject ranges
query_widths <- width(query[queryHits(hits)])
subject_widths <- width(subject[subjectHits(hits)])



# Filter overlaps with at least threshold overlap
keep <- (overlap_widths / query_widths >= threshold) | 
        (overlap_widths / subject_widths >= threshold)



# Subset the overlaps
threshold_hits <- hits[keep]


# Extract the filtered intersected ranges
threshold_overlap_ranges <- pintersect(query[queryHits(threshold_hits)], 
                                          subject[subjectHits(threshold_hits)])



# check the result
threshold_overlap_ranges

## GRanges object with 2 ranges and 2 metadata columns:
##     seqnames    ranges strand |     score       hit
##        <Rle> <IRanges>  <Rle> | <numeric> <logical>
##   a     chr1   150-200      + |        10      TRUE
##   b     chr1   300-400      + |        20      TRUE
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

TagMatrix

• TagMatrix : peaks在基因组区域内的位置分布，并为每个区域(bins)生成一个矩阵

  TagMatrix 的行代表 不同的peaks。

  TagMatrix 的列代表区域被划分的 bins

• 假设你对基因的 转录起始位点（TSS） 进行分析，并设定 upstream = 3000、downstream = 3000，且 nbin = 800：

  • 这 6000 bp 被划分为 800 个 bin，即 每个 bin 代表 6000/800 = 7.5 bp 的区间。

  • TagMatrix[i, j] 代表 第 i 个 peak 在第 j 个 bin 里的信号值（如读取数或加权值)

• getTagMatrix() can get TagMatrix from GRanges object

tagMatrix <- getTagMatrix(
  peak_list[[1]],       # peak file or GRanges object
  weightCol = "score",  # column name of peak as weight
  TxDb = TxDb,          # TxDb object
  upstream = 500,
  downstream = 500,
  verbose = TRUE,       # print message or not
  ignore_strand = FALSE,# 是否忽略链信息(默认区分正负链)
  
  nbin = 500,          # window region 被划分成多少个 bins（小区间）（用于计算信号密度）
                       # 值越大，分辨率越高，但可能会导致计算时间增加
  
  
  #windows,            # 接提供一个基因组区域(window region)的 GRanges 对象
                       # 适用于分析某个基因的外显子、内含子等特定window region的峰值分布
                       # 如果提供了 windows,那么 upstream 和 downstream 参数会被忽略）
  
  type = "start_site", # 选择绘制区域("start_site"：TSS, Transcription Start Site，
                       #               "end_site"：TES, Transcription End Site
                       #               "body" : Gene Body (整个基因体))
  
  by = "gene"         # 按'gene', 'transcript', 'exon', 'intron' , '3UTR' , '5UTR', 'UTR'进行分析
)

## >> binning method is used...2025-03-24 16:33:38
## >> preparing start_site regions by gene... 2025-03-24 16:33:38
## >> preparing tag matrix by binning...  2025-03-24 16:33:38

# Draw a table
kable(tagMatrix[c(1:10),c(1:10)], format = "html") %>%
  kable_styling(bootstrap_options = c("striped",     # 为表格添加斑马条纹（隔行变色）
                                      "hover",       # 在鼠标悬停时高亮当前行
                                      "condensed",   # 紧凑模式
                                      "responsive")) # 自适应布局

Heatmap and peak profiling

• you can use tagHeatmap() to create heatmap with TagMatrix

  tagHeatmap(tagMatrix)

  

• you can use peakHeatmap() to create heatmap with peak file or GRanges object

  peakHeatmap(
    peak_list[[1]],    # peak file or GRanges object
    weightCol = "score",  # column name of peak as weight
    TxDb = TxDb,        # TxDb object
    upstream = 500,     # upstream position  ,rel(0.2)
    downstream = 500,   # downstream position
    xlab = "",         # xlab
    ylab = "",         # ylab
    title = NULL,      
    palette = NULL,    # 颜色方案
    verbose = FALSE,   # print message or not
    by = "gene",       # 按'gene', 'transcript', 'exon', 'intron' , '3UTR' , '5UTR', 'UTR'进行分析
  
    type = "start_site", # 选择绘制区域("start_site"：TSS, Transcription Start Site，
                         #        "end_site"：TES, Transcription End Site
                         #        "body" : Gene Body (整个基因体))
  
    nbin = 500,          # window region 被划分成多少个 bins（小区间）（用于计算信号密度）
                         # 值越大，分辨率越高，但可能会导致计算时间增加
  
    ignore_strand = FALSE, # 是否忽略链信息（默认区分正负链）
  
    #windows,           # 接提供一个基因组区域(window region)的 GRanges 对象
                        # 适用于分析某个基因的外显子、内含子等特定window region的峰值分布
                        # 如果提供了 windows,那么 upstream 和 downstream 参数会被忽略）
  
  
    ncol = NULL,       # 指定输出的图像列数
    nrow = NULL        # 指定输出的图像行数
  )

  ## >> binning method is used...2025-03-24 16:34:04

  

• peakHeatmap_multiple_Sets() can compare of the peak heatmap over two regions

  #提取所有转录本的基因信息
  txdb1 <- transcripts(TxDb)
  
  #仅选择前 10,000 个基因的 5' UTR，避免计算量过大
  txdb2 <- unlist(fiveUTRsByTranscript(TxDb))[1:10000,]
  
  region_list <- list(geneX = txdb1, geneY = txdb2)
  
  peakHeatmap_multiple_Sets(peak_list[[1]],
                            upstream = 500,
                            downstream = 500,
                            by = c("geneX","geneY"),
                            type = "start_site",
                            TxDb = region_list,
                            nbin = 500,
                            verbose = FALSE)

  ## >> binning method is used...2025-03-24 16:34:27

  

• peak_Profile_Heatmap() can plot heatmap and peak profiling together

  peak_Profile_Heatmap(peak = peak_list[[1]],
                       upstream = 500,
                       downstream = 500,
                       by = "gene",
                       type = "start_site",
                       TxDb = TxDb,
                       nbin = 500,
                       verbose = FALSE)

  ## >> binning method is used...2025-03-24 16:35:07
  ## >> binning method is used...2025-03-24 16:35:11

  

• peak_Profile_Heatmap() is also support exploring several regions with heatmap and peak profiling

  #提取所有转录本的基因信息
  txdb1 <- transcripts(TxDb)
  
  #仅选择前 10,000 个基因的 5' UTR，避免计算量过大
  txdb2 <- unlist(fiveUTRsByTranscript(TxDb))[1:10000,]
  
  region_list <- list(geneX = txdb1, geneY = txdb2)
  
  
  peak_Profile_Heatmap(peak = peak_list[[1]],
                       upstream = 500,
                       downstream = 500,
                       by = c("geneX","geneY"),
                       type = "start_site",
                       TxDb = region_list,
                       nbin = 500,
                       verbose = FALSE)

  ## >> binning method is used...2025-03-24 16:35:33
  ## >> binning method is used...2025-03-24 16:35:39

  

Average Profile of peaks binding to TSS region

• plotAvgProf() can draw Average Profile of peaks with TagMatrix

• parametets:conf(confidence interval): 图中灰色的地方，conf = 0.95 表示 期望95%的样本结果落入这个区间

  plotAvgProf(tagMatrix, 
              xlim=c(1, 500),   # Need to be smaller than and match your nbin 
              xlab="Genomic Region (5'->3')", 
              ylab = "Read Count Frequency" ,
              conf = 0.95)

  ## >> plotting figure...             2025-03-24 16:36:04 
  ## >> Running bootstrapping for tag matrix...        2025-03-24 16:36:08

  

• plotAvgProf2() can draw Average Profile of peaks with peak file or GRanges object

  plotAvgProf2(
    peak_list[[1]],    # peak file or GRanges object
    weightCol = NULL,  # column name of peak as weight
    TxDb = TxDb,
    upstream = 500,
    downstream = 500,
    xlab = "Genomic Region (5'->3')",
    ylab = "Peak Count Frequency",
    conf = 0.95,        # confidence interval
    facet = "none",     # 分组："none"（默认）：不分组,"row"：按行分组,"column"：按列分组
    free_y = TRUE,      # 是否允许 Y 轴自由缩放  
    verbose = FALSE,    #  print message or not
    ignore_strand = FALSE,
  )

  ## >> plotting figure...             2025-03-24 16:36:09 
  ## >> Running bootstrapping for tag matrix...        2025-03-24 16:36:15

  

Pie plot

plotAnnoPie(peakAnno)

Bar plot

plotAnnoBar(peakAnno)

Ven pie plot

User may interested to view the full annotation with their overlap

vennpie(peakAnno)

Upsetplot()

We extend UpSetR to view full annotation overlap. User can user upsetplot() function.

upsetplot(peakAnno)

## Warning: Removed 8 rows containing non-finite outside the scale range
## (`stat_count()`).

We can combine vennpie with upsetplot() by setting vennpie = TRUE.

upsetplot(peakAnno, vennpie=TRUE)

## Warning: Removed 8 rows containing non-finite outside the scale range
## (`stat_count()`).

Visualize distribution of TF-binding loci relative to TSS

plotDistToTSS(peakAnno,
              title="Distribution of transcription factor-binding loci\nrelative to TSS")

• X 轴：代表结合位点的百分比（Binding sites (%)），即不同区间内的 TF 结合位点占比

Profile of several peak data binding to TSS region

• plotAvgProf(), tagHeatmap() and plotPeakProf() can accept a list of tagMatrix and visualize profile or heatmap among several Peaks

• plotAvgProf2(), peakHeatmap() and plotPeakProf2() can accept a list of bed files or GRange list and perform the same task

• the above function usage see [Heatmap and peak profiling] and [Average Profile of peaks binding to TSS region]

plotPeakProf2(peak_list, upstream = 3000, downstream = 3000, conf = 0.95,
              by = "gene", type = "start_site", TxDb = TxDb,
              facet = "row", nbin = 800)

## >> binning method is used...2025-03-24 16:36:25
## >> preparing start_site regions by gene... 2025-03-24 16:36:25
## >> preparing tag matrix by binning...  2025-03-24 16:36:25 
## >> binning method is used...2025-03-24 16:36:31
## >> preparing start_site regions by gene... 2025-03-24 16:36:31
## >> preparing tag matrix by binning...  2025-03-24 16:36:31 
## >> Running bootstrapping for tag matrix...        2025-03-24 16:36:45 
## >> Running bootstrapping for tag matrix...        2025-03-24 16:36:54

peak annotation comparision

• The plotAnnoBar() and plotDistToTSS() can also accept input of a named list of annotated peaks

# Use lapply to process each element
peakAnnoList <- lapply(peak_list, annotatePeak, TxDb=TxDb,
                       tssRegion=c(-3000, 3000), verbose=FALSE)

plotAnnoBar(peakAnnoList)

plotDistToTSS(peakAnnoList)

Functional profiles comparison

The clusterProfiler package provides compareCluster() function for comparing biological themes among gene clusters, and can be easily adopted to compare different peak experiments.

# Use lapply to process each element to get gene set
genes = lapply(peakAnnoList, function(i) as.data.frame(i)$geneId)

# enrichment analysis
compKEGG <- compareCluster(geneCluster   = genes,
                           fun           = "enrichKEGG",
                           pvalueCutoff  = 0.05,
                           pAdjustMethod = "BH",
                           organism = "rno",)

# dotplot
dotplot(compKEGG, showCategory = 15, title = "KEGG Pathway Enrichment Analysis")

Overlap of peaks and annotated genes

• vennplot accept a list of object, can be a list of GRanges or a list of vector.

• The detailed Overlap of peaks you should see [Peak overlap]

# a list of GRanges
vennplot(peak_list)

# a list of vector
genes= lapply(peakAnnoList, function(i) as.data.frame(i)$geneId)
vennplot(genes)

Shuffle genome coordination

• The shuffle() function to randomly permute the genomic locations of ChIP peaks defined in a genome which stored in TxDb object

• shuffle() 随机化 ChIP-seq 峰的基因组坐标，但确保它们仍然位于 基因组的有效区域。这个方法常用于计算 背景分布，进而评估 ChIP-seq 峰的 富集显著性

p <- GRanges(seqnames=c("chr1", "chr3"),
             ranges=IRanges(start=c(1, 100), end=c(50, 130)))

shuffle(p, TxDb=TxDb)

## GRanges object with 2 ranges and 0 metadata columns:
##       seqnames              ranges strand
##          <Rle>           <IRanges>  <Rle>
##   [1]     chr1   79002751-79002800      *
##   [2]     chr3 102133046-102133076      *
##   -------
##   seqinfo: 2 sequences from an unspecified genome; no seqlengths

Peak overlap enrichment analysis

enrichPeakOverlap() 工作原理

• 计算真实重叠峰值数

  计算 queryPeak和 targetPeak之间的实际重叠峰数量。

• 生成随机背景分布

  使用 shuffle() 函数 随机化 queryPeak 的基因组坐标，生成 nShuffle 组 随机峰数据。

  对每组随机峰集，与 targetPeak 计算重叠峰值数，从而构建 背景分布（null distribution）。

• 计算富集显著性

  统计实际峰值重叠数 是否远高于随机背景，用 p 值 评估其显著性。

  采用 Benjamini-Hochberg（BH）方法进行 多重假设检验校正来校正 p 值。

Overlap_analysis=enrichPeakOverlap(
  queryPeak=peak_list[[1]],       # bed file or GRanges object
  targetPeak=peak_list,  # target bed file(s) or folder or a list of GRanges objects
  TxDb= TxDb,           # TxDb objects
  pAdjustMethod="BH",   # pvalue adjustment method
  nShuffle=50,          # shuffle numbers
  chainFile=NULL,       # chain file for liftOver(基因轴坐标转换)
  pool=FALSE,           # whether Merge targetPeak (multiple ChIP seq peaks)
  mc.cores = detectCores()-1,  # number of cores
  verbose= FALSE)  # whether Print detailed information       

• Error: The input of list of GRanges for targetPeak here will result in different target peak p-values being the same. If we input bed file (s) or folder , this problem will not occur.

• To speed up the compilation of this vignettes, we only set nShuffle to 50 as an example for only demonstration. User should set the number to 1000 or above for more robust result

• Parameter chainFile are chain file name for mapping the targetPeak to the genome version consistent with queryPeak when their genome version are different.

  • This creat the possibility of comparison among different genome version and cross species

  • The detailed chainFile you should see [Chain file]

kable(Overlap_analysis, format = "html") %>%
  kable_styling(bootstrap_options = c("striped",#为表格添加斑马条纹（隔行变色）
                                      "hover",       # 在鼠标悬停时高亮当前行
                                      "condensed",   # 紧凑模式
                                      "responsive")) # 自适应布局 

• qSample is the name of queryPeak

• qLen is the the number of peaks in queryPeak

• N_OL is the number of overlap between queryPeak and targetPeak

GEO data collection

There are many ChIP seq data sets that have been published and deposited in GEO database

ChIPseekercollect about 17,000 bed files deposited in GEO

• User can use getGEOspecies() to get a summary based on speices.

  GEOspecies <-getGEOspecies()

  # Draw a table
  reactable(GEOspecies, 
            striped = TRUE,    # 斑马条纹
            highlight = TRUE,  # 鼠标悬停高亮
            bordered = TRUE,   # 添加表格边框
            filterable = TRUE, # 允许按列筛选数据 
            searchable = TRUE, # 允许全局搜索表格内容
            resizable = TRUE   # 允许拖拽调整列宽
  )

• User can use getGEOgenomeVersion() to get a summary based on genome version

  GEOgenomeVersion <- getGEOgenomeVersion()

  # Draw a table
  reactable(GEOgenomeVersion, 
            striped = TRUE,    # 斑马条纹
            highlight = TRUE,  # 鼠标悬停高亮
            bordered = TRUE,   # 添加表格边框
            filterable = TRUE, # 允许按列筛选数据 
            searchable = TRUE, # 允许全局搜索表格内容
            resizable = TRUE   # 允许拖拽调整列宽
  )

• User can access the detail information by getGEOInfo(), for each genome version.

  rn5 <- getGEOInfo(genome="rn5", simplify=TRUE)

  If simplify is set to FALSE

  extra information including source_name, extract_protocol, description, data_processing and submission_date.

  reactable(rn5, 
            striped = TRUE,    # 斑马条纹
            highlight = TRUE,  # 鼠标悬停高亮
            bordered = TRUE,   # 添加表格边框
            filterable = TRUE, # 允许按列筛选数据 
            searchable = TRUE, # 允许全局搜索表格内容
            resizable = TRUE,  # 允许拖拽调整列宽
            height = 800
  )

Download GEO ChIP data sets

ChIPseeker provide function downloadGEObedFiles() to download all the bed files of a particular genome

# creat rn5 folder
dir.create("rn5", showWarnings = FALSE, recursive = TRUE)  

这里作者提供的downloadGEObedFiles() 代码有问题，因为GEO 服务器可能已经不支持 FTP 访问，建议改用 HTTPS，所以我对代码进行了修改，创建了个新函数downloadGEObedFilesHTTPS()

downloadGEObedFilesHTTPS <- function(genome,destDir) {
    # Obtaining genomic information
    info <- getGEOInfo(genome)
    
    # Get supplementary_file column
    fnames <- as.character(info$supplementary_file)
    
    # Replace FTP with HTTPS
    fnames_https <- gsub("^ftp://", "https://", fnames)
    
    # Generate target file path
    destfiles <- file.path(destDir, basename(fnames_https))
    names(destfiles) <- NULL
    
    # Traverse and download files
    for (i in seq_along(fnames_https)) {
        
      # 检查 destfiles[i]是否已存在，文件不存在，则执行下载
      if (!file.exists(destfiles[i])) {
            tryCatch(
                download.file(fnames_https[i], 
                              destfile = destfiles[i], 
                              mode = "wb"),#以二进制模式下载，确保文件不会因换行符问题损坏。
                error = function(e) message(fnames_https[i], ": file not found and skip")
            )
        }
    }
}

Overlap significant testing

• After download the bed files from GEO, we can pass them to for testing the significant of overlap

• Signifcant overlap can be use to generate hypothesis of cooperative regulation

• We can use enrichPeakOverlap() function to assess Statistical testing . detail see [Statistical testing of ChIP seq overlap]

• In different genome versions, parameter chainFile that can pass a chain file and liftOver the targetPeak to the genome version consistent with queryPeak